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We develop a model to study the thermal expansion of surfaces, wherein phonon frequencies 
are obtained from ab initio total energy calculations. Anharmonic effects are treated exactly 
in the direction normal to the surface, and within a quasiharmonic approximation in the plane 
of the surface. We apply this model to the Ag(lll) and Al(lll) surfaces, and find that our 
calculations reproduce the experimental observation of a large and anomalous increase in the 
surface thermal expansion of Ag(lll) at high temperatures^. Surprisingly, we find that this 
increase can be attributed to a rapid softening of the in-plane phonon frequencies, rather than 
due to the anharmonicity of the out-of-plane surface phonon modes. This provides evidence for a 
new mechanism for the enhancement of surface anharmonicity. A comparison with Al(lll) shows 
that the two surfaces behave quite differently, with no evidence for such anomalous behavior on 
Al(lll). 
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The equilibrium lattice constant of a crystal is determined by a balance between the 
various attractive and repulsive forces present within the solid. When the crystal is cleaved 
to form surfaces, this balance is destroyed, and the atoms at the surface therefore relax 
either inwards (which is the usual situation for most metal surfaces) or outwards. 

However, upon heating the crystal, these relaxations may change dramatically as a func- 
tion of temperature. The phenomenon of bulk thermal expansion is, of course, a familiar 
one: upon heating the crystal, the lattice constant changes, reflecting the anharmonicity of 
the interatomic potentials. Similarly, asymmetries and anharmonicities in surface phonon 
vibrations will be manifested in a change in interlayer displacements as a function of tem- 
perature T. Due to the different asymmetries present at the surface and in the bulk, it is 
possible that these two quantities (bulk lattice constant and interlayer distances near the 
surface) may change quite differently upon increasing T. 

Indeed, it has long been realized that one would expect all such measures of anharmonic- 
ity (e.g., coefficients of thermal expansion a, mean squared displacements (MSDs) of atoms, 
and the rate of change of phonon frequency with temperature) to be larger at surfaces than 
in bulk crystals^ for two reasons: (i) the breaking of symmetry due to the presence of the 
surface makes the interlayer potential more asymmetric at the surface than in the bulk, 
increasing the size of the odd terms in a Taylor series expansion of the energy in powers of 
atomic displacements; (ii) atomic displacements may be greater at the surface than in the 
bulk, thus increasing the relative magnitudes of the higher-order anharmonic terms in this 
series expansion. 

Interestingly, studies on metal surfaces have shown that such enhancements of surface 
anharmonicity, when observed, are strongly dependent on both the element and the ori- 
entation of the surface. Early experiments and calculations^^ on metal surfaces suggested 
that measures of anharmonicity (such as the coefficient of thermal expansion) are typically 
at most three times larger at the surface than in the bulk, and as a result do not affect 
surface properties drastically. However, recent experiments have shown a huge enhancement 
on a few surfaces: for Ni(OOl)! and Pb(110)@ at high temperatures, the surface coefficient 
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of thermal expansion as, defined by: 

as = (di2)-\dd 12 /dT), (f) 

(where dyi is the interlayer spacing between the first two planes of atoms at the surface) is 
10 to 20 times larger than the bulk coefficient of thermal expansion Ob, while for Cu(110) 
the surface MSDs are up to six times greater than the bulk msdH. 

Perhaps the most interesting case is that of Ag(lll): recent experiments^ show that the 
contraction of dyi is dramatically reversed for the Ag(lll) surface as the crystal is heated: 
Up to T « 670 K, d\2 is indeed contracted by ~ 2.5% relative to the bulk separation ds] 
but upon increasing T further, <i 12 increases much more rapidly than <i B does, so that by 
1150 K it is expanded by ~ 10%. Correspondingly, as becomes more than ten times as 
large as aj. Such a large effect is especially unexpected for a close-packed metal surface, 
and contradicts the conventional expectation that more open surfaces should exhibit larger 
surface anharmonicity. 

Lewis@ has carried out EAM simulations to investigate the thermal behavior of Ag(lll). 
The results of these simulations differ significantly from those reported experimentally: the 
surface layer relaxes inwards at all temperatures, and as is less than twice as large as 

Is this disagreement between experiment and calculation on Ag(lll) due to inadequacies 
of the EAM potentials? Or could it be a sign of some hitherto undetected surface phase 
transition? 

To study these questions, we have investigated the harmonic and anharmonic proper- 
ties of Ag(lll) and Al(lll) by performing ah initio calculations using density functional 
theory. Fully separable norm-conserving pseudopotentials@ were used in our calculations, 
together with a plane wave basis set with an energy cut-off of 60 Ry (20 Ry for Al), and the 
local-density approximation with Ceperley-Alder exchange-correlation^. We note that the 
relatively high cut-off of the plane-wave basis set is necessary in order to obtain a good de- 
scription of anharmonic effects, even though a lower cut-off may suffice to describe harmonic 
properties adequately. 



Before performing surface calculations, we first verified that the harmonic and anhar- 
monic properties of the bulk materials are satisfactorily described by these pseudopotentials. 
Some of these results, such as the optimal values of the lattice constant and the bulk modu- 
lus, are presented in Table 1. This table also contains the calculated values of two measures 
of anharmonicity: the pressure derivative of the bulk modulus, and the Griineisen parame- 
ters 7, which describe how the phonon frequencies vary upon changing the lattice constant. 
Note that the experimental value of 7 is an average over all bulk modes, whereas the theo- 
retical values were calculated separately for each band, and at a sample wave-vector along 
the [111] direction, between the zone-center T and the zone-edge L. (These particular bulk 
modes were chosen because they project on to the zone-center of the surface Brillouin zone, 
and can therefore be regarded as analogous to the surface zone-center vibrations that we 
will later investigate for the (111) surface.) 

The surface calculations were performed using a repeated slab geometry consisting of six 
atomic layers separated by a vacuum layer of the same thickness. The k-point sets used to 
sample reciprocal space consisted of a uniform grid centered on the T point and containing 
seven points in the irreducible part of the Brillouin zone for the undistorted surface; the 
number of k-points was correspondingly increased upon breaking symmetries by distorting 
the lattice in order to calculate phonon frequencies. Convergence of calculated anharmonic 
quantities with respect to energy cut-off, number of k-points and number of layers was 
carefully tested for. 

Our strategy is to compute static energies and phonon frequencies by performing self- 
consistent calculations at T = K, and then extend our results to finite temperatures by 
using a quasiharmonic approximation. 

In order to obtain a qualitative understanding of the mechanisms in operation, and a 
first estimate of the size of surface thermal effects, we use a simple model of the lattice 
dynamics of the surface, considering only three phonon modes, in all of which the topmost 
layer of the slab moves as a whole - i.e., we assume that the displacements are confined to 
the first layer of atoms at the surface, and consider only those modes with zero wave-vector. 
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We consider one mode in which the surface atoms vibrate normal to the surface plane (along 
the z direction), and two modes (along x = [110] and y = [112]) in which they vibrate in 
the plane of the surface. ( We emphasize that these displacements do not correspond to 
any of the actual normal modes of vibration of the surface slab, but may be considered 
as indicative of the strengths of the various force-constants that would appear in the true 
dynamical matrix of the system. This proviso should be kept in mind when we refer to 
"modes" and "phonons" in the rest of this paper.) 

We first computed the change in the total energy of the Ag(lll) slab upon varying 
the first two interlayer separations d 12 and 0^23, and maintaining the fee stacking of the 
bulk crystal. This not only provides the static interlayer potential, but is also equivalent 
to simulating the vibrational mode along z. Our result for the dependence on <i 12 of the 
first interlayer potential is plotted in Fig. 1; it is clearly asymmetric about the minimum at 
d\2 = 2.30 A. We found that allowing for the relaxation of c/23 does not have a significant 
impact on the results for the close-packed (111) surfaceB. For the results presented in Fig. 1 
and the rest of this paper, 0^23 is therefore fixed at the bulk interlayer separation of 2.34 A. 

To see how this anharmonicity of the interlayer potential is manifested at finite tempera- 
tures, we consider a one-dimensional quantum oscillator vibrating in the interlayer potential 
of Fig. 1. A numerical solution of the Schrodinger equation for this problem furnishes the 
eigenstates and eigenvalues of such an oscillator, and the mean displacement (di2) n i n the 
n-th eigenstate is obtained by computing the expectation value of the displacement operator 
in each state. The average value at a finite temperature T is then obtained by weighting 
these results with the corresponding partition function. 

Our results for c/i 2 (T) obtained from this procedure (see the open circles in Fig. 3) 
indicate a modest enhancement in as relative to of ~ 1.7, which is much smaller than 
that measured experimentally. 

We next performed frozen-phonon calculations to study the behavior of the two in-plane 
modes in our model. At each value of <ii 2 , the atoms in the surface layer were displaced 
along first the x and then the y direction, and the total energy was computed for a series 

5 



of displacements up to ±0.15 A. The curvature of the resulting plots of energy versus dis- 
placement gives the mode frequency. We find that the frequency of these in-plane modes 
decreases surprisingly rapidly upon increasing c? 12 ; these results are plotted in Fig. 2. 

In analogy to the usual bulk Griineisen parameter 7#, we can define a surface Griineisen 
parameter: 

ls = -dHcu)/dln(d 12 ). (2) 

The value of the 75 extracted from our results for the in-plane modes (between 5.5 and 7.5) is 
significantly larger than the corresponding value of of 2.39 that we obtain for a bulk vibration 
with an analogous pattern of displacements. Thus, the frequently made assumption^ that 7 
is approximately equal for surface and bulk modes is clearly invalid in this case. 

The surface can thus reduce its vibrational free energy significantly by expanding out- 
wards, though such an outward expansion would be accompanied by an increase in the static 
energy. The optimal value of dyi is determined by minimizing the free energy function!: 

F(d 12 ,T) = E stat (d l2 ) + J2Kiddi2,T), (3) 

i 

at each temperature T . Here, -E l s tat(^i2) is the static interlayer potential plotted in Fig. 1, and 
F* ih (di2,T) is the vibrational free energy corresponding to vibrations in the i-th direction, 
which, in the quasiharmonic approximation, is given byi: 

FL(di2,T) = k B Tln{2sinh{^^-)}- (4) 

where ks and h are Boltzmann's constant and Planck's constant respectively. The frequency 
of the mode polarized along the i-th direction, when the first interlayer separation is fixed 
at di2, is denoted by u;j(d 12 ). The sum in Eq. (3) runs over all the bands of the phonon 
spectrum (averaged over the entire surface Brillouin zone); in our case we approximate it 
by a sum over the 3 zone-center vibrational patterns that we have considered, which are 
polarized along the x,y and z directions respectively. The variation of ^ (c? 12 ) is obtained 
directly from our frozen-phonon calculations for the two in-plane modes (see Fig. 2). For the 
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out-of-plane mode, we compute F* ih (di2,T) numerically, in such a way as to reproduce the 
exact result for di 2 (T) that we have already obtained by solving the Schrodinger equation 
when only the mode along z is present. 

Our final result for d\ 2 (T) in the presence of all three modes, obtained by minimizing 
the free energy expression given by Eq. (3), is shown by the filled circles in Fig. 3 (note 
that the temperatures are normalized with respect to the melting temperature T m ). The 
experimentally measured data pointsEI are also plotted; both the experimental and theoretical 
curves display the same features: there is little or no change up to about T/T m = 0.5, i.e., 
T = 617 K, beyond which there is a rapidly increasing trend towards outwards relaxation of 
the surface layer. At T/T m = 0.85, i.e., T = 1049 K, we find that the surface layer is relaxed 
outwards by about 15%, whereas the experiments show an outwards relaxation of ~ 7.5%. 
Given the simplicity of our model, and the large experimental error bars, this is as good an 
agreement as we can hope for. 

The increasing slope of d 12 (T) reflects a flattening in the minimum of the free-energy 
curve, and the rapid increase in di 2 (T) at high T is a precursor to the development of a 
saddle-point instability in the free-energy curve, similar to that which has been obtained in 
studies of the surface melting of copper surfaces^. We emphasize that the rapid decrease of 
^i{di 2 ) for the in-plane modes is crucial for obtaining the large outwards expansion; if, for 
example, this rate of decrease were to be halved, the maximum outwards expansion would 
be drastically reduced to about 2%. In the limit of high T (when all modes are excited), the 
value of di 2 (T) is no longer sensitive to the absolute scale of u>i, but is instead controlled by 
7s, which is a normalized indicator of how rapidly u>i falls off with increasing d± 2 . 

To check whether such behavior is universal or a peculiar property of Ag(lll), we re- 
peated the same calculations on bulk aluminum and Al(lll). Our results for the relaxation 
of du for Al(lll) are also plotted in Fig. 3, and it is obvious that there is no evidence for 
a dramatically increased surface expansion on Al(lll). We have also performed calcula- 
tions on Cu(lll)@ which show that the behavior of Cu(lll) is intermediate between that 
of Ag(lll) and Al(lll), which has been confirmed by very recent experiments^. 



Our results indicate that in addition to the two well known sources of enhanced surface 
anharmonicity that we have already mentioned, a third (and hitherto neglected) effect is 
more important in causing the dramatic enhancement in surface anharmonicity on Ag(lll): 
Not only do interlayer potentials at the surface tail off rapidly with increasing z , they 
simultaneously become much flatter in the xy plane - in other words, the operative effect is 
not so much a decrease in the absolute magnitudes of interlayer potentials at the surface, but 
a reduction in their corrugation parallel to surface. As a consequence, those surface phonon 
modes in which atomic dispacements have significant components in the surface plane soften 
rapidly upon increasing interlayer separations; this drives the outermost layer of atoms to 
expand outwards at high temperatures. 

An accurate description of the in-plane corrugation clearly requires taking into account 
the correct distribution of the electronic charge density (and the resulting chemical bond 
formation) at the surface. It seems plausible that when rebonding effects become significant, 
the EAM (which essentially ignores the relaxation of atomic charge densities and the rehy- 
bridization of electronic states) may fail; this may be why Lewis! did not observe a large 
outwards expansion in his simulations. 

The rapid decrease in the corrugation of the interlayer potential is evident in Fig. 4, 
where we have plotted the differences in energies when the outermost layer of atoms occupies 
various stacking sites. Note that (%) the flattening occurs more rapidly for Ag(lll) than 
Al (111) (ii) upon allowing for the lighter mass of Al atoms, the effective corrugation relevant 
for phonon frequencies is actually larger for Al(lll) than for Ag(lll). Both these factors 
contribute to the enhancement in 7 of the in-plane top-layer modes and thus the larger 
thermal expansion of Ag(lll). 

Such a decrease in the corrugation of the substrate potential would also tend to favor a 
top-layer reconstruction of the type that has been observed on Au(lll)@or Pt(lll)@, where 
the substrate potential is too weak to prevent a densification of atoms in the topmost layer. 
However, experiments apparently show no evidence of such a reconstruction on Ag(lll)S; 
further calculations of surface stresses and the strength of intralayer couplings should help 
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clarify the situation. 

Our conclusion that the enhancement in surface expansion arises mainly from in-plane 
vibrations is supported by the experimental observation that the amplitude of in-plane 
vibrations on Ag(lll) rises faster than the magnitude of out-of-plane vibrations@. Exper- 
iments on other surfaces, e.g., Cu(OOl), have detected in-plane vibrational amplitudes that 
are larger than out-of-plane amplitudesil, and it is interesting to speculate whether this 
(counter-intuitive) result arises from the same cause, i.e., from a rapid softening of in-plane 
frequencies due to thermal expansion. 

Of course, the calculation we have presented above does contain certain approximations: 
we have restricted ourselves to surface vibrations at the zone-center, we approximate the 
true electronic exchange-correlation functional by the LDA, and we use a quasiharmonic 
form for vibrational free energy. We should also allow for expansion in the plane of the 
surface, but this should in fact reinforce the softening of the in-plane modes. The validity 
of some of these approximations can be tested by performing ab initio molecular dynamics 
simulations in which all anharmonic contributions are fully included automatically, and in 
which no assumptions are made about the polarizations of bulk and surface vibrations; work 
in this direction is in progress. 

We expect that our numerical results will change slightly upon including other surface 
phonon modes and allowing for dispersion through the surface Brillouin zone. Surface vi- 
brations at the zone-center are sensitive only to the strength of interlayer force-constants 
coupling atoms at the surface to atoms in layers below; however, the frequency of a phonon 
at arbitrary wave- vector also depends on the intralayer force-constants, and one may expect 
these to be less sensitive to di 2 ; we may therefore have over-estimated the degree of surface 
anharmonicity by restricting ourselves to the zone-center. However, we note that measure- 
ments of surface phonon frequency- shifts suggest that the degree of anharmonicity remains 
approximately constant through a surface phonon band0. 

In conclusion, we have developed a simple model for the study of surface thermal expan- 
sion. By applying it to Ag(lll), we have demonstrated that the anomalously large surface 



thermal expansion of Ag(lll) can be attributed to a rapid softening of in-plane vibrational 
modes (related to a rapid flattening of the corrugation of the interlayer potential) upon 
increasing interlayer distances. We have shown that a similar scenario does not, however, 
lead to significant enhancement on Al(lll). 
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Table 1 



Quantity 


Ag 


Al 


do (A) 


4.06 (4.09) 


3.94 (4.05) 


S(MBar) 


i.22 ^.w; 


0.80 f0.75j 


B' 


5.94 f5.£2j 


4.22 (5.36) 


7(T) 


2.39 (246) 


2.12 




2.95 


2.05 ^gj«; 



Table 1: Calculated values of bulk properties of Ag and Al: a = the equilibrium lattice 
constant at T = 0, 5=bulk modulus, B' = pressure derivative of B, j(T) = Griineisen 
parameter for transverse bulk mode along [111], with wave- vector 2/3 of the way between 
the zone-center and zone-edge, ^y(L) = Griineisen parameter for longitudinal bulk mode at 
the same wave- vector. Experimental values at room temperature, obtained from Reference^, 
are given in parentheses. 
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Fig.l: Static interlayer potential between the two outermost layers of Ag(lll). The hori- 
zontal lines and open circles indicate the energy eigenvalues E n and the mean displacements 
{d\2)n respectively of an Ag atom vibrating in this potential. 
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Fig. 2: Energy TiuJi of the top-layer in-plane modes of Ag(lll), as a function of the first 
interlayer separation d\2- Circles and squares indicate modes polarized along the x- and y- 
directions respectively. 
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Fig. 3: Contraction/expansion of d 12 relative to ds, as a function of normalized tem- 
perature T/T m - our calculations for Ag(lll) with out-of-plane mode only (open circles), 
and all three modes (filled circles); experiment 1 on Ag(lll) (solid line with error bars); our 
calculation for Al(lll) with all three modes (stars). T rn is the bulk melting temperature 
{T^f = 1234 K, T£ = 933 K). 
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Fig. 4: Decrease in the corrugation of the interlayer potential with increasing interlayer 
separation rf 12 : the graphs show the increase in surface energy when the outermost layer of 
atoms occupies an atop site or bridge site instead of the favored fee hollow site, for Ag(lll) 
(filled circles) and Al(lll) (stars). Arrows indicate the equilibrium value (neglecting zero- 
point vibrations) of di 2 at T = K. 
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